{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Часть 1. Описание набора данных и признаков. Цель исследования\n",
    "#### Основной набор данных:\n",
    "\n",
    "Этот набор данных содержит цены продажи жилья для округа Кинг, с административным центром в г. Сиэтл за период с мая 2014 года по май 2015 года.\n",
    "Краткая информация из Википедии об округе Кинг (англ. King County)  \n",
    "Крупные города: Сиэтл (окружной центр и самый большой город), Белвью, Такома\n",
    "Население округа: 2 188 649 (на 2017 год)\n",
    "Территория:\t\n",
    " - Общая площадь\t2,307 кв. миль (5,975 km2)\n",
    " - Площадь земель\t2,116 кв. миль (5,480 km2)\n",
    " - Площадь водоемов\t191 кв. миль (495 km2)\n",
    "Сиэтл (англ. Seattle) — крупнейший город на северо-западе США и в штате Вашингтон, крупный морской порт. Расположен между системой заливов Пьюджет и озером Вашингтон.\n",
    "\n",
    "#### Информация о наборе данных\n",
    "Адресс нахождения данных: https://www.kaggle.com/harlfoxem/housesalesprediction/data\n",
    "\n",
    "<b>id</b> - Номер обеъкта\n",
    "\n",
    "\n",
    "<b>date</b> - Дата продажи объекта неджвижимости\n",
    "\n",
    "\n",
    "<b>price</b>- Цена продажи\n",
    "\n",
    "\n",
    "<b>bedrooms</b> - Кол-во спален\n",
    "\n",
    "\n",
    "<b>bathrooms</b> - Кол-во ванных+сан.узлов\n",
    "\n",
    "\n",
    "<b>sqft_living</b>- Жилая площадь дома\n",
    "\n",
    "\n",
    "<b>sqft_lot</b> - Общая площадь участка\n",
    "\n",
    "\n",
    "<b>'floors'</b> - Кол-во этажей в доме\n",
    "\n",
    "\n",
    "<b>\"waterfront\"</b> - Вид на водоем (Да/Нет)\n",
    "\n",
    "\n",
    "<b>view</b> - Был просмотрен (не очень понял что тут имелось ввиду)\n",
    "\n",
    "\n",
    "<b>condition</b> - Насколько хорошо состояние (в целом)\n",
    "\n",
    "\n",
    "<b>grade</b> - Общий класс(грейд), присвоенный обеъкту неджвижимости, основанный на системе классификации Кинг Кантри\n",
    "\n",
    "\n",
    "<b>sqft_above</b> - Площадь объекта недвижимости, не включая подвал(если есть)\n",
    "\n",
    "\n",
    "<b>sqft_basement</b> - Площадь подвала (если есть)\n",
    "\n",
    "<b>yr_built</b> - Год постройки\n",
    "\n",
    "\n",
    "<b>yr_renovated</b> - Год реновации (если была)\n",
    "\n",
    "\n",
    "<b>zipcode</b> - Почтовый индекс\n",
    "\n",
    "\n",
    "<b>lat</b> - Широта\n",
    "\n",
    "\n",
    "<b>long</b> - Долгота\n",
    "\n",
    "\n",
    "<b>sqft_living15</b> - Жилая площадь в 2015 году (подразумевает некоторые ремонтные работы/ перепланировку). Это может повлиять или не повлиять на общую площадь территории\n",
    "\n",
    "\n",
    "<b>sqft_lot15</b> - Площадь общей территории в 2015 году(подразумевает некоторые ремонтные работы/ перепланировку)\n",
    "\n",
    "\n",
    "#### Доплнительные данные взяты отсюда:\n",
    "\n",
    "https://moto.data.socrata.com/dataset/King-County-Sheriff-s-Office/4h35-4mtu - Информация о преступлениях в округе Кинг\n",
    "\n",
    "https://data.kingcounty.gov - Информационный портал округа Кинг(несколько разных отчетов по тематикам, Школы,Клиники,Магазины и общепит из которых я вытаскивал местонаходжение объектов и сами объекты)\n",
    "\n",
    "Данные в агрегированном виде и основной датасет: https://cloud.mail.ru/public/DTDE/zax4jHwSL\n",
    "\n",
    "#### Цель исследования: \n",
    "1. Исследовать данные по продажам недвижимости округа Кинг.\n",
    "\n",
    "2. Построить модель наиболее точно предсказывающую стоимость объектов недвижимости на основании имеющихся данных и за счет сбора дополнительных.\n",
    "\n",
    "#### Первичный план исследования:\n",
    "1. Для ислледования и поиска лучшей модели я будут использовать на первом этапе регрессию Lasso и RandomForestRegressor с дефолтными или почти дефолтными настройками.\n",
    "2. Метрики качества R2  и RMSLE - буду смотреть на обе, так удобнее. Первая характеризует некую абстрактную аккуратность, вторая превязана к реальным деньгам.\n",
    "3. Запуск моделей буду проводить после каждого серьезного шага улучшений и добавлений новых признаков, чтобы следить за тем, что общий тренд исследования стремится к улучшению.\n",
    "4. В конце попробую градиентный бустинг и далее блендинг и стеккинг."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Часть 2. Первичный анализ данных\n",
    "#### Импорт библиотек"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Импорт всех нужных библиотек\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "\n",
    "from sklearn.model_selection import train_test_split, GridSearchCV, cross_val_predict, cross_val_score\n",
    "from sklearn.base import BaseEstimator, ClassifierMixin\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score\n",
    "from sklearn.linear_model import Lasso, Ridge\n",
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.cluster import KMeans\n",
    "import xgboost as xgb\n",
    "import lightgbm as lgb\n",
    "\n",
    "from scipy import stats\n",
    "from scipy.special import inv_boxcox, boxcox1p\n",
    "from scipy.stats.mstats import gmean\n",
    "from scipy.stats import norm, skew , boxcox\n",
    "from scipy.sparse import csr_matrix\n",
    "\n",
    "pd.set_option('display.max_columns', 100)\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Чтение данных"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('kc_house_data.csv', parse_dates  =['date'])\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Данные не имеют пропусков, что уже хорошо."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Сразу отложу данные для валидации исследования"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ind = len(data)*0.8\n",
    "work_data = data.ix[:ind]\n",
    "validation_data = data.ix[ind:]\n",
    "data.shape, work_data.shape, validation_data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Зададим первичный Baseline, который будем потом улучшать\n",
    "Перед анализом данных, поскольку сам дата сет без пропусков и ошибок, я построю простую модель RandomForest и использую его метрику качества score - по умолчанию R2 - коэф. детерминации, которую он считает на отложенной выборке. Дополнительно буду смотреть на среднеквадратичную логарифмическую ошибку (достаточно частая метрика в пердсказании стоимости недвижимости). Также посмотрю как они (R2 и RSMLE) зависимы друг от друга практически. RandomForest также выбираю из практических соображений, чтобы сразу посмотреть на значимость признаков"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Функция среднеквадратичной логарифмической ошибки\n",
    "def rmsle(y_true,y_pred):\n",
    "    assert len(y_true) == len(y_pred)\n",
    "    return np.square(np.log(y_pred + 1) - np.log(y_true + 1)).mean() ** 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train = work_data.drop([ 'id','price', 'date'], axis=1)\n",
    "target = work_data['price']\n",
    "X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)\n",
    "\n",
    "\n",
    "clf = RandomForestRegressor(random_state=17)\n",
    "clf.fit(X_train.values, y_train)\n",
    "\n",
    "y_pred = clf.predict(X_test.values)\n",
    "clf.score(X_test.values, y_test), rmsle(y_test, y_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Сразу посмотрим какие признаки алгоритм оранжировал по степени важности\n",
    "feature_importance = pd.DataFrame()\n",
    "feature_importance['feature'] = df_train.columns\n",
    "feature_importance['importances'] = clf.feature_importances_\n",
    "feature_importance.sort_values('importances', ascending=False).head(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Часть 3. Визуальный анализ данных\n",
    "\n",
    " Посмотрим на график зависимости цены от первых 5 самых важных признаков из базовго алгоритма, а также\n",
    " их зависимость между собой"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set()\n",
    "cols = feature_importance['feature'][feature_importance['importances']> 0.031]\n",
    "sns.pairplot(work_data[cols], size = 2.5)\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Из графиков видно, что на высокий грейд линейно влияет общая жилая площадь, а также немного жилая площадь по оценке на 2015 год (учитывающая реновации)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Посмотрим как распределены продажи по zipcode\n",
    "\n",
    "sort_data = work_data.groupby('zipcode', sort = 'id').count().reset_index()\n",
    "plt.figure(figsize=(15,10))\n",
    "ax= sns.barplot(x=sort_data['zipcode'], y=sort_data['bedrooms'])\n",
    "plt.xticks(rotation= 90)\n",
    "plt.xlabel('Zipcode')\n",
    "plt.ylabel('Кол-во объектов')\n",
    "plt.title(\"Распределение кол-ва проданных объектов относительно их районов, выделенных по zipcode\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В некоторых районах недвижимости продается в разы больше чем в других. Невысокиие продажи могут говорить как правило о двух противоположных тенеденциях. Дорогой район и там состав жильцов не меняется. Бедный район и там мало кто покупает. Посмотрим на других срезах, подтвердится ли гипотеза."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Посмотрим как распределена средняя стоимость проданных объектов относительно их районов, выделенных по zipcode\n",
    "\n",
    "sort_data_pz = work_data.groupby('zipcode', sort = 'id').mean().reset_index()\n",
    "plt.figure(figsize=(15,10))\n",
    "ax= sns.barplot(x=sort_data_pz['zipcode'], y=sort_data_pz['price'])\n",
    "plt.xticks(rotation= 90)\n",
    "plt.xlabel('Zipcode')\n",
    "plt.ylabel('Средняя стоимость')\n",
    "plt.title(\"Распределение средней стоимости проданных объектов относительно их районов, выделенных по zipcode\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Предположение:\n",
    "Видно что есть несколько индексов(4), среденяя стоимость жилья возле которых превышает 1 млн. долларов. По всей видимости это какие-то элитные районы. Можно создать дополнительно 3 признака класса жилья по стоимости. \n",
    "Частичное подтверждение гипотезы по индексу 98039 - продается мало объектом и средняя стоимостьсвыше 2млн. долларов. - богатый район.\n",
    "\n",
    "98039 - Медина, Вашингтон, США - входит в топ 10 самых дорогих районов в США, ср. стоимость недвижимости (по данным из интернета за 2015 год) около 2 998 000 долларов, в нашем графике около 2 100 000, но у нас все что реально продалось, а аналитические журналы часто аппелируют средней заявленной стартовой ценой продавца."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sort_data_yr = work_data.groupby('yr_built').mean().reset_index()\n",
    "plt.figure(figsize=(15,10))\n",
    "ax= sns.barplot(x=sort_data_yr['yr_built'], y=sort_data_yr['price'])\n",
    "plt.xticks(rotation= 90)\n",
    "plt.xlabel('Год постройки')\n",
    "plt.ylabel('Средняя стоимость')\n",
    "plt.title(\"Распределение средней стоимости относительно года постройки\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Интересно, что на графике видно, что средняя стоимость постоек до 1941 года и после 1994 выше, чем в период между этими датами. Это может свидетельствовать о таких факторах, как:\n",
    "- Старые исторические здания в популярных районах (возможно исторические центры) ценятся выше.\n",
    "- Новые постройки обычно строятся в экономически интересных для потребителей районах, а как следствиее спрос диктует цену\n",
    "    \n",
    "Т.е. возможно, постройки до 41 и после 94 - это строительство примерно в одних и тех же районах.\n",
    "Нужно построить график распределения по годам построек по индексам."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15,10))\n",
    "sort_data_zy = work_data.groupby(['zipcode','yr_built']).count().reset_index()\n",
    "sns.boxplot(x=\"zipcode\", y=\"yr_built\", data=sort_data_zy)\n",
    "plt.xticks(rotation= 90)\n",
    "plt.xlabel('Индекс')\n",
    "plt.ylabel('Год постройки')\n",
    "plt.title(\"Распределение возраста объектов недвижимости в засисимости от индекса\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Предположение\n",
    "Очевидно, что есть более \"молодые районы\" и более \"старинные\", также можно предположить, что индексы с \"широким ящиком\" это районы близкие к историческому центру городов или они и являются центрами, также из теории урбанизации можно сделать вывод что такие районы ближе к воде (к историческому порту), если город находится у судоходного водоема.\n",
    "\n",
    "Проверил индекс 98102 - это почти центр г. Сиэтл и близкий к воде - т.е. предположение подтвердилось"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Выделим из даты месяц продажи и посмотрим на распределение продаж по месяцам\n",
    "work_data['month_sold'] = work_data['date'].apply(lambda ts: ts.month)\n",
    "validation_data['month_sold'] = validation_data['date'].apply(lambda ts: ts.month)\n",
    "\n",
    "plt.figure(figsize=(10,7))\n",
    "sns.countplot(x=work_data.month_sold)\n",
    "plt.xlabel('Месяц продажи')\n",
    "plt.ylabel('Кол-во объектов')\n",
    "plt.title('График кол-ва продаж по месяцам',color = 'blue',fontsize=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видна определенная сезоннасть весной продают больше, чем зимой и т.д."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Посмотрим на расперделение по грейду\n",
    "plt.figure(figsize=(10,7))\n",
    "sns.countplot(x=work_data.grade)\n",
    "plt.xlabel('Грейд')\n",
    "plt.ylabel('Кол-во объектов')\n",
    "plt.title('График распределения грейда среди всех проданных объектов',color = 'blue',fontsize=15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Построим график зависимости цены на недвижимость в зависимости от самых значимых признаков (значимость больше 0.01)\n",
    "\n",
    "for i in feature_importance['feature'][feature_importance['importances']> 0.01]:\n",
    "    fig, ax = plt.subplots(figsize=(10, 10))\n",
    "    ax.scatter(x = work_data[i], y = work_data['price'])\n",
    "    plt.ylabel('Price', fontsize=13)\n",
    "    plt.xlabel(i, fontsize=13)\n",
    "    plt.show()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Из полезного на графике №2 наблюдается 1 выброс, который стоит удалить, дабы он не повлиял на модели, чувствительные к выбросам, если мы в итоге такие будем использовать\n",
    "\n",
    "По графикам зависимости от широты и долготы видно, что в дата сете присуствует порядка 10 объектов, у которых местоположение сильно влисяет на стоимость, возможно это условная Luxury недвижимость в элитных районах, где продается не так много объектов ежегодно."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Удалим выброс и построим график для проверки снова\n",
    "work_data = work_data[work_data['sqft_living'] < 13000]\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "ax.scatter(x = work_data['sqft_living'], y = work_data['price'])\n",
    "plt.ylabel('Price', fontsize=13)\n",
    "plt.xlabel('Sqft living Area', fontsize=13)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Полученные наблюдения распределений буду использовать для создания новых пороговых признаков"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Посмотрим на количество уникальных значений признаков\n",
    "for col in work_data.select_dtypes(exclude=['datetime64']).columns:\n",
    "    print([col], len(work_data[col].unique()))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Из интересного здесь.\n",
    "\n",
    "Категориальные признаки:\n",
    " - zipcode , floors, yr_built, yr_renovated, bedrooms, bathrooms, waterfront\n",
    " \n",
    "Ранговые:\n",
    " - grade, condition\n",
    " \n",
    "эту информацию можно будет использовать для генерации новых признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Построим стандартную корреляционную матрицу для наглядности по всем числовым признакам"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corrmat = work_data.corr()\n",
    "f, ax = plt.subplots(figsize=(15, 10))\n",
    "sns.heatmap(corrmat, square=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Из интересного здесь можно увидеть несколько интересных корреляционных групп.\n",
    "    \n",
    "Достаточно сильная корреляция между собой и ценой у:\n",
    "    \n",
    "    - bedroom, bathroom, sqft_living\n",
    "    \n",
    "    - waterfront, view\n",
    "    \n",
    "    - grade, sqft_above"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Теперь посмотрим что из себя вообще представляет ценнобразование - наша целевая метрика, с точки зрения распределения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(8, 8)\n",
    "\n",
    "sns.distplot(work_data['price'] , fit=norm);\n",
    "(mu, sigma) = norm.fit(data['price'])\n",
    "print( '\\n mu = {:.2f} and sigma = {:.2f}\\n'.format(mu, sigma))\n",
    "\n",
    "#Построим график распределения\n",
    "plt.legend(['Нормальное распределение ($\\mu=$ {:.2f} and $\\sigma=$ {:.2f} )'.format(mu, sigma)],\n",
    "            loc='best')\n",
    "plt.ylabel('Частота цены')\n",
    "plt.title('Распределение цены')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Дополнительно построим QQ-plot\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "res = stats.probplot(work_data['price'], plot=plt)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "На 1-ом графике видно, что распределение цены имеет правый уклон, что на самомм деле достаточно типично, по моему опытуЭ для рынка недвижимости США. По крайней мере 5-6 датасетов, что мне встречались имели похожий тренд.\n",
    "\n",
    "Уже из этих грфиков, в принципе можно предположить,что линейные модели, которые заточены под нормальное распределение, будут иметь, скорее всего не очень высокий результат без преобразования данных."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Проверим гипотезу"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train = work_data.drop([ 'id','price', 'date'], axis=1)\n",
    "target = work_data['price']\n",
    "X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)\n",
    "scalar = StandardScaler()\n",
    "clf = Lasso(random_state=17)\n",
    "lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])\n",
    "lasso_base.set_params(Lasso__alpha=1).fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = lasso_base.predict(X_test)\n",
    "print(\" R2_score: {}, RMSLE: {}\".format(lasso_base.score(X_test, y_test),rmsle(y_test, y_pred)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Унылость результата подтверждает гипотезу. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Посчитаем значение лямбда для price, чтобы преобразовать значение признака и построим график снова"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sk_date = work_data['price'] - min(work_data['price']) + 1e-5\n",
    "_, lmbda_price  = boxcox(sk_date)\n",
    "target = boxcox1p(work_data['price'],  lmbda_price)\n",
    "#Также проведем туже операцию для данных на валидации\n",
    "target_val = boxcox1p(validation_data['price'],  lmbda_price)\n",
    "print (lmbda_price)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(8, 8)\n",
    "sns.distplot(target , fit=norm);\n",
    "(mu, sigma) = norm.fit(work_data['price'])\n",
    "print( '\\n mu = {:.2f} and sigma = {:.2f}\\n'.format(mu, sigma))\n",
    "\n",
    "#Построим график распределения\n",
    "plt.legend(['Нормальное распределение ($\\mu=$ {:.2f} and $\\sigma=$ {:.2f} )'.format(mu, sigma)],\n",
    "            loc='best')\n",
    "plt.ylabel('Частота цены')\n",
    "plt.title('Распределение цены')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Дополнительно построим QQ-plot\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "res = stats.probplot(target, plot=plt)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Почти прекрасно"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### На правах \"Эксперимента\" - любопытства ради"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Попробуем нормализавать распределение через геометрическое среднее\n",
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(8, 8)\n",
    "geom = (work_data['price']**0.104 - 1 )/ 0.104 * gmean(work_data['price'], axis=0)**(0.104 - 1)\n",
    "sns.distplot(geom , fit=norm);\n",
    "\n",
    "plt.ylabel('Частота цены')\n",
    "plt.title('Распределение цены')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В данном случае, просто проверил формулу, которую видел в учебниках по статистике"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Посмотрим на распределение других количественных признаков с широким диапазоном значений\n",
    "columns_disc = ['sqft_living', 'sqft_above', 'sqft_basement', 'sqft_living15', 'sqft_lot15', 'sqft_lot']\n",
    "for col in columns_disc:\n",
    "    fig, ax = plt.subplots()\n",
    "    fig.set_size_inches(8, 8)\n",
    "    sns.distplot(work_data[col] , fit=norm);\n",
    "    (mu, sigma) = norm.fit(work_data[col])\n",
    "    \n",
    "\n",
    "    #Построим график распределения\n",
    "    plt.legend(['Нормальное распределение ($\\mu=$ {:.2f} and $\\sigma=$ {:.2f} )'.format(mu, sigma)],\n",
    "            loc='best')\n",
    "    plt.ylabel('Частота')\n",
    "    plt.title('Распределение')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Очевидно что все вышеуказанные колличественные признаки также как и цена имеют правый уклон.\n",
    "\n",
    "_lot15, _lot, _basement - имеют много нулевых значений. Возможн стоит дополнительно добавить преобразованные дискретные признаки \"Да/Нет\" от них в набор данных"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Часть 4. Преобразование исходных признаков."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "skewed_feats = work_data[columns_disc].apply(lambda x: skew(x)).sort_values(ascending=False)\n",
    "#col = list(data.drop(['date', 'id', 'zipcode','price', 'lat', 'long'], axis=1).columns)\n",
    "#skewed_feats = data[col].apply(lambda x: skew(x)).sort_values(ascending=False)\n",
    "print(\"\\nУклон количественных признаков: \\n\")\n",
    "skewness = pd.DataFrame({'Skew' :skewed_feats})\n",
    "skewness.head(6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для нормально распределенных признаков значение skew должно быть равно 0 или быть близко к нему."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для того, чтобы избежать во входных данных появления отрицательных или равных нулю значений, всегда будем находить минимальное значение входной последовательности и вычитать его из каждого ее элемента, дополнительно осуществляя сдвиг на небольшую величину, равную 1e-5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Вычислим оптимальное lambda для каждого\n",
    "skewed_features = skewness.index\n",
    "lambd = []\n",
    "for feat in skewed_features:\n",
    "    work_data[feat] = work_data[feat] - min(work_data[feat]) + 1e-5\n",
    "    _, lmbda = boxcox(work_data[feat])\n",
    "    lambd.append(lmbda)\n",
    "    print(feat, lmbda)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum([abs(x) for x in lambd])/6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "##### Возьму итоговую общую лямбду =  0.12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "lam = 0.12\n",
    "for feat in skewed_features:\n",
    "    work_data[feat] = work_data[feat] - min(work_data[feat]) + 1e-5\n",
    "    work_data[feat] = boxcox1p(work_data[feat], lam)\n",
    "    validation_data[feat] = validation_data[feat] - min(validation_data[feat]) + 1e-5\n",
    "    validation_data[feat] = boxcox1p(validation_data[feat], lam)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in columns_disc:\n",
    "    fig, ax = plt.subplots()\n",
    "    fig.set_size_inches(8, 8)\n",
    "    sns.distplot(work_data[col] , fit=norm);\n",
    "    (mu, sigma) = norm.fit(work_data[col])\n",
    "    \n",
    "\n",
    "    #Построим график распределения\n",
    "    plt.legend(['Нормальное распределение ($\\mu=$ {:.2f} and $\\sigma=$ {:.2f} )'.format(mu, sigma)],\n",
    "            loc='best')\n",
    "    plt.ylabel('Частота')\n",
    "    plt.title('Распределение')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Кроме площади цокольного этажа почти все признаки стали более нормально распределенными. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Повторное контрольное обучение модели на испраленных признаках"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train = work_data.drop(['id','price', 'date'], axis=1)\n",
    "X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)\n",
    "\n",
    "\n",
    "clf = RandomForestRegressor(random_state=17)\n",
    "rnd_base = Pipeline([('scalar', scalar), ('rnd', clf)])\n",
    "rnd_base.fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = rnd_base.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(rnd_base.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Есть улучшение относительно исходных данных (0.8504728683309748, 0.19446848336861014) - улучшились оба показателя. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим что произошло с важностью признаков."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_importance = pd.DataFrame()\n",
    "feature_importance['feature'] = df_train.columns\n",
    "feature_importance['importances'] = clf.feature_importances_\n",
    "feature_importance.sort_values('importances', ascending=False).head(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Интересно. Широта обогнала жилую площадь."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Посмотрим как улучшилась регрессионная модель.\n",
    "df_train = work_data.drop(['id','price', 'date'], axis=1)\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)\n",
    "scalar = StandardScaler()\n",
    "clf = Lasso(random_state=17)\n",
    "lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])\n",
    "lasso_base.set_params(Lasso__alpha=1).fit(X_train.values, y_train)\n",
    "                     \n",
    "y_pred = lasso_base.predict(X_test.values)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(lasso_base.score(X_test, y_test), \n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Показатели регрессионной модели ухудшились.\n",
    "\n",
    "Было: R2_score: 0.711325562859866, RMSLE: 0.39307021980748347\n",
    "\n",
    "Вообще это странно - такое ухудшение. Возможно, модель просто недообучается с таким высоким значением alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Проверим коэффициенты\n",
    "print(\"Коэффициенты: {}\\n всего коэффициентов: {}\\n всего коэффициентов отличныхот 0: {}\".format(clf.coef_, len(clf.coef_),\n",
    "                                                                                                 len(clf.coef_[clf.coef_ != 0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Что и следовало доказать, высокий коэффициент регуляризации занижает много весов, оставляя только 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Уменьшим регуляризацию\n",
    "lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])\n",
    "lasso_base.set_params(Lasso__alpha=0.1).fit(X_train.values, y_train)\n",
    "                     \n",
    "y_pred = lasso_base.predict(X_test.values)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(lasso_base.score(X_test, y_test), \n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))\n",
    "print(\"Коэффициенты: {}\\n всего коэффициентов: {}\\n всего коэффициентов отличныхот 0: {}\".format(clf.coef_, len(clf.coef_),\n",
    "                                                                                                 len(clf.coef_[clf.coef_ != 0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Построим кривую кросс-валидации\n",
    "cv_scores, holdout_scores = [], []\n",
    "alphas = np.linspace(0.001, 1.0, 15)\n",
    "\n",
    "for i in alphas:\n",
    "\n",
    "    lasso = Lasso(alpha=i)\n",
    "    cv_scores.append(np.mean(cross_val_score(lasso, X_train, y_train, cv=3)))\n",
    "    lasso.fit(X_train, y_train)\n",
    "    holdout_scores.append(r2_score(y_test, lasso.predict(X_test)))\n",
    "\n",
    "plt.figure(figsize=(8, 8))\n",
    "plt.plot(alphas, cv_scores, label='CV')\n",
    "plt.plot(alphas, holdout_scores, label='holdout')\n",
    "plt.title('Easy task. Lasso fails')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Лучше оставить такой alpha = 0.1, так гораздо больше признаков работает."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Часть 5. Создание новых признаков из начальных данных\n",
    "Займемся поиском точек роста качества модели за счет инжиниринга признаков.\n",
    "В части 3, где мы проводили визуальный анализ было обнаружено несколько интересныз закономерностей, которые мы инкрементируем в модель данных и посмотрим, как они повлияют на качество обучения."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Добавим дополнительные дискретные признаки для колличественных, где большая часть имеет нулевые значения\n",
    "work_data['_lot'] = 0\n",
    "work_data['_lot15'] = 0\n",
    "work_data['_basement'] = 0\n",
    "work_data['_lot'][work_data['sqft_lot']>0] = 1\n",
    "work_data['_lot15'][work_data['sqft_lot15']>0] = 1 \n",
    "work_data['_basement'][work_data['sqft_basement']>0] = 1\n",
    "\n",
    "# Выделю 3 группы объектов исходя из их стоимости\n",
    "work_data['chip'] = 0\n",
    "work_data['mid'] = 0 \n",
    "work_data['exp'] = 0\n",
    "work_data['chip'][work_data['price']<= 375000] = 1\n",
    "work_data['mid'][(work_data['price']> 375000)&(work_data['price']<= 1000000)] = 1 \n",
    "work_data['exp'][work_data['price']> 1000000] = 1\n",
    "\n",
    "#Выделю сезоны продажи недвижимости\n",
    "work_data['sum'] = 0\n",
    "work_data['spr'] = 0\n",
    "work_data['aut'] = 0\n",
    "work_data['win'] = 0\n",
    "work_data['sum'][(work_data['month_sold']>= 6) & (work_data['month_sold']<= 8)] = 1\n",
    "work_data['spr'][(work_data['month_sold']>= 3) & (work_data['month_sold']<= 5)] = 1\n",
    "work_data['aut'][(work_data['month_sold']>= 9) & (work_data['month_sold']<= 11)] = 1\n",
    "work_data['win'][(work_data['month_sold']>= 1) & (work_data['month_sold']<= 2) | (work_data['month_sold']== 12)] = 1\n",
    "\n",
    "# Выделю клас жилья через грейд\n",
    "work_data['low'] = 0\n",
    "work_data['average'] = 0 \n",
    "work_data['high'] = 0\n",
    "work_data['low'][work_data['grade']<= 5] = 1\n",
    "work_data['average'][(work_data['price']> 5)&(work_data['price']<= 9)] = 1 \n",
    "work_data['high'][work_data['price']> 9] = 1\n",
    "\n",
    "#Выделю 3 группы исходя из года постройки\n",
    "work_data['old'] = 0\n",
    "work_data['old_to_hihg'] = 0 \n",
    "work_data['young'] = 0\n",
    "work_data['old'][work_data['yr_built']<= 1942] = 1\n",
    "work_data['old_to_hihg'][(work_data['yr_built']> 1942)&(work_data['yr_built']<= 1994)] = 1 \n",
    "work_data['young'][work_data['yr_built']> 1994] = 1\n",
    "\n",
    "#Добавим 4 новых признака, 2 характеризующих  \"возраст\" дома, наличие подвала \"Да\\Нет\", была ли реноваци \"Да\\нет\"\n",
    "work_data['years_old'] = 2015 - work_data['yr_built']\n",
    "work_data['second_young'] = 2015 - work_data['yr_renovated']\n",
    "work_data['second_young'][work_data['second_young'] == 2015] = 0\n",
    "work_data['basement'] = 0\n",
    "work_data['basement'][work_data['sqft_basement'] >0] = 1\n",
    "work_data['was_renov'] = 0\n",
    "work_data['was_renov'][work_data['yr_renovated'] >0] = 1\n",
    "#Напишем простую функцию для лучшего разделения признаков, характеризующих размеры\n",
    "def quadratic(feature):\n",
    "    work_data[feature+'2'] = work_data[feature]**2\n",
    "    \n",
    "quadratic('sqft_living')\n",
    "quadratic('sqft_living15')\n",
    "quadratic('sqft_lot')\n",
    "quadratic('sqft_above')\n",
    "quadratic('sqft_basement')\n",
    "\n",
    "# Проделю все тоже самое и для контрольной выборки\n",
    "validation_data['_lot'] = 0\n",
    "validation_data['_lot15'] = 0\n",
    "validation_data['_basement'] = 0\n",
    "validation_data['_lot'][validation_data['sqft_lot']>0] = 1\n",
    "validation_data['_lot15'][validation_data['sqft_lot15']>0] = 1 \n",
    "validation_data['_basement'][validation_data['sqft_basement']>0] = 1 \n",
    "\n",
    "validation_data['chip'] = 0\n",
    "validation_data['mid'] = 0 \n",
    "validation_data['exp'] = 0\n",
    "validation_data['chip'][validation_data['price']<= 375000] = 1\n",
    "validation_data['mid'][(validation_data['price']> 375000)&(validation_data['price']<= 1000000)] = 1 \n",
    "validation_data['exp'][validation_data['price']> 1000000] = 1\n",
    "\n",
    "validation_data['sum'] = 0\n",
    "validation_data['spr'] = 0\n",
    "validation_data['aut'] = 0\n",
    "validation_data['win'] = 0\n",
    "validation_data['sum'][(validation_data['month_sold']>= 6) & (validation_data['month_sold']<= 8)] = 1\n",
    "validation_data['spr'][(validation_data['month_sold']>= 3) & (validation_data['month_sold']<= 5)] = 1\n",
    "validation_data['aut'][(validation_data['month_sold']>= 9) & (validation_data['month_sold']<= 11)] = 1\n",
    "validation_data['win'][(validation_data['month_sold']>= 1) & \\\n",
    "                       (validation_data['month_sold']<= 2) | (validation_data['month_sold']== 12)] = 1\n",
    "\n",
    "validation_data['low'] = 0\n",
    "validation_data['average'] = 0 \n",
    "validation_data['high'] = 0\n",
    "validation_data['low'][validation_data['grade']<= 5] = 1\n",
    "validation_data['average'][(validation_data['price']> 5)&(validation_data['price']<= 9)] = 1 \n",
    "validation_data['high'][validation_data['price']> 9] = 1\n",
    "\n",
    "validation_data['old'] = 0\n",
    "validation_data['old_to_hihg'] = 0 \n",
    "validation_data['young'] = 0\n",
    "validation_data['old'][validation_data['yr_built']<= 1942] = 1\n",
    "validation_data['old_to_hihg'][(validation_data['yr_built']> 1942)&(validation_data['yr_built']<= 1994)] = 1 \n",
    "validation_data['young'][validation_data['yr_built']> 1994] = 1\n",
    "\n",
    "# Тоже самое для валидации\n",
    "validation_data['years_old'] = 2015 - validation_data['yr_built']\n",
    "validation_data['second_young'] = 2015 - validation_data['yr_renovated']\n",
    "validation_data['second_young'][validation_data['second_young'] == 2015] = 0\n",
    "validation_data['basement'] = 0\n",
    "validation_data['basement'][validation_data['sqft_basement'] >0] = 1\n",
    "validation_data['was_renov'] = 0\n",
    "validation_data['was_renov'][validation_data['yr_renovated'] >0] = 1\n",
    "\n",
    "def quadratic(feature):\n",
    "    validation_data[feature+'2'] = validation_data[feature]**2\n",
    "    \n",
    "quadratic('sqft_living')\n",
    "quadratic('sqft_living15')\n",
    "quadratic('sqft_lot')\n",
    "quadratic('sqft_above')\n",
    "quadratic('sqft_basement')\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Сверим, не упустил ли какие признаки для валидации"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "work_data.shape, validation_data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Проверим прогресс"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train = work_data.drop(['id','price', 'date'], axis=1)\n",
    "X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)\n",
    "\n",
    "\n",
    "clf = RandomForestRegressor(random_state=17)\n",
    "rnd_base = Pipeline([('scalar', scalar), ('rnd', clf)])\n",
    "rnd_base.fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = rnd_base.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(rnd_base.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Есть неплохой прирост (R2_score: 0.8757049690628793, RMSLE: 0.18806603652440654)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = Lasso(random_state=17)\n",
    "lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])\n",
    "lasso_base.set_params(Lasso__alpha=0.1).fit(X_train.values, y_train)\n",
    "                     \n",
    "y_pred = lasso_base.predict(X_test.values)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(lasso_base.score(X_test, y_test), \n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Регрессия тоже растет"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Часть 6. Создание новых признаков за счет загрузки доп. региональных данных*. \n",
    "*Данные были подготовлены отдельно через обработку сырых таблиц и отчетности.\n",
    "\n",
    "Уровень преступности - агрегированан из отчета с сайта шерифа округа, где фиксируются все преступления по времени, по типу и месту\n",
    "\n",
    "Школы - агрегирован, если правильно помню, из отчета по вакцинации учебных заведений округа (нач. школы и дет. сады).\n",
    "\n",
    "Клиники - таблица с дислокацией мед. учереждений\n",
    "\n",
    "Магазины, кафе, рестораны - агрегирован из отчета испекторов санитарных служб, которые занимаюца лецензированием таких заведений.\n",
    "\n",
    "Сами данные брались из публичных источников и обрабатывались для выявления местонахождения объектов. Идея состоит в том, что близость к определенным объектам, наличие объектов для базовых потребностей(школа, больница, магазин), скопление объектов или событий(криминал) - все это влияет на конечную стоимость жилья.\n",
    "Проверим..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Уровень преступности \n",
    "crime = pd.read_csv('crimeKC.csv', index_col=0)\n",
    "# Местоположения по индексу образовательных учереждений (нач. школа и дет. сады)\n",
    "school = pd.read_csv('schoolKC.csv', index_col=0)\n",
    "# Местоположения по индексу медицинских учереждений\n",
    "clinic = pd.read_csv('clinicsKC.csv', index_col=0)\n",
    "# Местоположения по индексу продуктовых магазинов, кафе, ресторанов\n",
    "shop_horeca = pd.read_csv('shop_cafeKC.csv', index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crime.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "school.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clinic.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shop_horeca.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Новые данные представляют из себя 4 разреженных матрицы, которые я сгенерировал из различных публичных государственных источников округа Кинг"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Часть 6.1 Начнем с уровня преступности\n",
    "В данной таблице собраны суммарные происшествия за несколько лет по типам и по zipcode.  Несмотря на то что данные агрегированны за период больший периоду продаж недвижимости, они скорее отражают тренд, выделяя более \"преступные\" локации."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "work_data = work_data.merge(crime, how='left', on='zipcode')\n",
    "#Некоторые данные не будут иметь совпадений из таблицы crime\n",
    "work_data = work_data.fillna(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "work_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Проверим регрессию\n",
    "df_train = work_data.drop(['id', 'price', 'date'], axis=1)\n",
    "X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)\n",
    "\n",
    "scalar = StandardScaler()\n",
    "clf = Lasso(random_state=17)\n",
    "lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])\n",
    "lasso_base.set_params(Lasso__alpha=0.1).fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = lasso_base.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(lasso_base.score(X_test, y_test), \n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Напомню до этого результат был: R2_score: 0.859201120304148, RMSLE: 0.1985696581535856 - т.е. новые данные улучшили модель."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Сравним с RandomForest\n",
    "\n",
    "clf = RandomForestRegressor(random_state=17)\n",
    "rnd_base = Pipeline([('scalar', scalar), ('rnd', clf)])\n",
    "rnd_base.fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = rnd_base.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(rnd_base.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "RandomForestRegressor немного ухудшился от прошлого значения: R2_score: 0.9169665913882967, RMSLE: 0.15419440365697273"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Часть 6.2 Добавим клиники и больницы"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "work_data = work_data.merge(clinic, how='left', on='zipcode')\n",
    "work_data = work_data.fillna(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "df_train = work_data.drop(['id', 'price', 'date'], axis=1)\n",
    "X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)\n",
    "\n",
    "scalar = StandardScaler()\n",
    "clf = Lasso(random_state=17)\n",
    "lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])\n",
    "lasso_base.set_params(Lasso__alpha=0.1).fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = lasso_base.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(lasso_base.score(X_test, y_test), \n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "df_train = work_data.drop(['id','price', 'date'], axis=1)\n",
    "X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)\n",
    "\n",
    "\n",
    "clf = RandomForestRegressor(random_state=17)\n",
    "rnd_base = Pipeline([('scalar', scalar), ('rnd', clf)])\n",
    "rnd_base.fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = rnd_base.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(rnd_base.score(X_test, y_test), \n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Снова небольшой рост на обоих моделях."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Часть 6.3 Добавим школы и детские сады."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "work_data = work_data.merge(school, how='left', on ='zipcode')\n",
    "work_data = work_data.fillna(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "df_train = work_data.drop(['id', 'price', 'date'], axis=1)\n",
    "X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)\n",
    "\n",
    "scalar = StandardScaler()\n",
    "clf = Lasso(random_state=17)\n",
    "lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])\n",
    "lasso_base.set_params(Lasso__alpha=0.1).fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = lasso_base.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(lasso_base.score(X_test, y_test), \n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "clf = RandomForestRegressor(random_state=17)\n",
    "rnd_base = Pipeline([('scalar', scalar), ('rnd', clf)])\n",
    "rnd_base.fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = rnd_base.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(rnd_base.score(X_test, y_test), \n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Часть 6.4 Добавим школы и детские сады."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "work_data = work_data.merge(shop_horeca, how='left', on ='zipcode')\n",
    "work_data = work_data.fillna(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "df_train = work_data.drop(['id', 'price', 'date'], axis=1)\n",
    "X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)\n",
    "\n",
    "scalar = StandardScaler()\n",
    "clf = Lasso(random_state=17)\n",
    "lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])\n",
    "lasso_base.set_params(Lasso__alpha=0.1).fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = lasso_base.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(lasso_base.score(X_test, y_test), \n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "clf = RandomForestRegressor(random_state=17)\n",
    "rnd_base = Pipeline([('scalar', scalar), ('rnd', clf)])\n",
    "rnd_base.fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = rnd_base.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(rnd_base.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Предварительный вывод:\n",
    "\n",
    "1. Работа с признаками позволила улучшить качество базовых алгоритмов.\n",
    "2. Случайный лес практически из коробки существенно лучше и быстрее справлялся с задачей на сырых данных.\n",
    "3. Регрессия намного чувствительнее к нормальности распределения, как таргета, так и признаков и с дефолтовой регуляризацией работает для данного набора данных плохо.\n",
    "4. Если скрость работы алгоритма будет критичным аспектом, то Lasso - предпочтительнее, при условии получения приемлемого качества за счет настройки гиперпараметров( я в нее пока верю:)).\n",
    "5. Необходимо проверить алгоритмя с оптимальными гиперпараметрами, а аткже проверить и другие алгоритмы для данной задачи.\n",
    "6. Проверка моделей после каждого нового набора данных показала, что все новые данные улучшали качество модели, что в свою очередь подтвердило мою изначальную гипотезу о влиянии данных факторов на стоимость объектов недвижимости."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Часть 7. Выбор модели. Кросс-валидация и настройка гиперпараметров модели."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Воспользуемся функцией одной из лекции курса\n",
    "def plot_with_err(x, data, **kwargs): \n",
    "    mu, std = data.mean(1), data.std(1)\n",
    "    lines = plt.plot(x, mu, '-', **kwargs)\n",
    "    plt.fill_between(x, mu - std, mu + std, edgecolor='none',\n",
    "    facecolor=lines[0].get_color(), alpha=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "steps = [('scalar', StandardScaler()),\n",
    "         ('Lasso', Lasso(random_state=17))]\n",
    "pipeline = Pipeline(steps)\n",
    "parameters = {'Lasso__alpha': [0.001, 0.01, 0.1, 1]}\n",
    "\n",
    "cvLasso = GridSearchCV(pipeline,param_grid=parameters,cv=3)\n",
    "cvLasso.fit(X_train.values,y_train)\n",
    "print(\"Accuracy: {}\".format(cvLasso.score(X_test, y_test)))\n",
    "print(\"Tuned Model Parameters: {}\".format(cvLasso.best_params_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvLasso.cv_results_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Посмотрим, что у нас получилось на кросс-валидации для Lasso.\n",
    "plt.figure(figsize=(8, 8))\n",
    "test_results = np.array([cvLasso.cv_results_['split0_test_score'],cvLasso.cv_results_['split1_test_score'],\n",
    "                     cvLasso.cv_results_['split2_test_score']]).T;\n",
    "train_results = np.array([cvLasso.cv_results_['split0_train_score'],cvLasso.cv_results_['split1_train_score'],\n",
    "                      cvLasso.cv_results_['split2_train_score']]).T;\n",
    "plot_with_err([0.001, 0.01, 0.1, 1], test_results, label='validation scores')\n",
    "plot_with_err([0.001, 0.01, 0.1, 1], train_results, label='train scores')\n",
    "plt.xlabel('alphas'); plt.ylabel('R2')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Качество падает линейно с ростом коэффициента регуляризации. Регрессия очень не стабильна."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "steps = [('scalar', StandardScaler()),\n",
    "         ('RandomForestRegressor', RandomForestRegressor(random_state=17))]\n",
    "pipeline = Pipeline(steps)\n",
    "parameters = {'RandomForestRegressor__n_estimators': [50, 100, 150],\n",
    "             'RandomForestRegressor__max_depth' : [15, 20, 25]}\n",
    "\n",
    "cvRf = GridSearchCV(pipeline,param_grid=parameters,cv=3)\n",
    "cvRf.fit(X_train,y_train)\n",
    "print(\"Accuracy: {}\".format(cvRf.score(X_test, y_test)))\n",
    "print(\"Tuned Model Parameters: {}\".format(cvRf.best_params_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Внимание: при  n_estimators более 250 (однажды выставил такое значение) - Считалось около 3.5 часов. \n",
    "Случайный лес с повышением при увеличение числа деревьев и максимальной глубины, наращивает качество, но становится очень медленным и начинает проигрывать по скорости регрессии очень сильно. Вообще я запускал перебор только 2 раза, а признаки добавлял новые гораздо чаще, так что при запуске score будет отличным от текущего."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvRf.cv_results_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Посмотрим, что у нас получилось на кросс-валидации для RandomForest.\n",
    "plt.figure(figsize=(8, 8))\n",
    "test_results = np.array([cvRf.cv_results_['split0_test_score'],cvRf.cv_results_['split1_test_score'],\n",
    "                     cvRf.cv_results_['split2_test_score']]).T;\n",
    "train_results = np.array([cvRf.cv_results_['split0_train_score'],cvRf.cv_results_['split1_train_score'],\n",
    "                      cvRf.cv_results_['split2_train_score']]).T;\n",
    "plot_with_err(['1-15_50', '2-15_100', '3-15_150', '4-20_50', '5-20_100', '6-20_150',\n",
    "               '7-25_50', '8-25_100', '9-25_150'], test_results, label='validation scores')\n",
    "plot_with_err(['1-15_50', '2-15_100', '3-15_150', '4-20_50', '5-20_100', '6-20_150',\n",
    "               '7-25_50', '8-25_100', '9-25_150'], train_results, label='train scores')\n",
    "plt.xlabel('depth-n_estiators'); plt.ylabel('R2')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RandomForest намного стабильнее Lasso, причем уже при n_estimators = 100 и глубине = 15 уже дает отличный результат"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "steps = [('scalar', StandardScaler()),\n",
    "         ('lgb.LGBMRegressor', lgb.LGBMRegressor(random_state=17))]\n",
    "pipeline = Pipeline(steps)\n",
    "parameters = {'lgb.LGBMRegressor__n_estimators': [700, 800,900],\n",
    "             'lgb.LGBMRegressor__max_depth' : [3, 4, 5],\n",
    "             'lgb.LGBMRegressor__learning_rate' : [0.03, 0.06, 0.09]}\n",
    "\n",
    "cvLgb = GridSearchCV(pipeline,param_grid=parameters,cv=3)\n",
    "cvLgb.fit(X_train,y_train)\n",
    "print(\"Accuracy: {}\".format(cvLgb.score(X_test, y_test)))\n",
    "print(\"Tuned Model Parameters: {}\".format(cvLgb.best_params_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvLgb.cv_results_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Посмотрим, что у нас получилось на кросс-валидации для lgb.\n",
    "plt.figure(figsize=(15, 15))\n",
    "test_results = np.array([cvLgb.cv_results_['split0_test_score'],cvLgb.cv_results_['split1_test_score'],\n",
    "                     cvLgb.cv_results_['split2_test_score']]).T;\n",
    "train_results = np.array([cvLgb.cv_results_['split0_train_score'],cvLgb.cv_results_['split1_train_score'],\n",
    "                      cvLgb.cv_results_['split2_train_score']]).T;\n",
    "plot_with_err(['11-0.03_3_700', '11-0.03_3_800', '11-0.03_3_900', '12-0.03_4_700', '12-0.03_4_800', '12-0.03_4_900',\n",
    "               '13-0.03_5_700', '13-0.03_5_800', '13-0.03_5_900', '21-0.06_3_700', '21-0.06_3_800', '21-0.06_3_900', \n",
    "               '22-0.06_4_700', '22-0.06_4_800', '22-0.06_4_900','23-0.06_5_700', '23-0.06_5_800', '23-0.06_5_900',\n",
    "               '31-0.09_3_700', '31-0.09_3_800', '31-0.09_3_900', '32-0.09_4_700', '32-0.09_4_800', '32-0.09_4_900',\n",
    "               '33-0.09_5_700', '33-0.09_5_800', '33-0.09_5_900'], test_results, label='validation scores')\n",
    "plot_with_err(['11-0.03_3_700', '11-0.03_3_800', '11-0.03_3_900', '12-0.03_4_700', '12-0.03_4_800', '12-0.03_4_900',\n",
    "               '13-0.03_5_700', '13-0.03_5_800', '13-0.03_5_900', '21-0.06_3_700', '21-0.06_3_800', '21-0.06_3_900', \n",
    "               '22-0.06_4_700', '22-0.06_4_800', '22-0.06_4_900','23-0.06_5_700', '23-0.06_5_800', '23-0.06_5_900',\n",
    "               '31-0.09_3_700', '31-0.09_3_800', '31-0.09_3_900', '32-0.09_4_700', '32-0.09_4_800', '32-0.09_4_900',\n",
    "               '33-0.09_5_700', '33-0.09_5_800', '33-0.09_5_900'], train_results, label='train scores')\n",
    "plt.xlabel('learning rate - depth - n_estiators'); plt.ylabel('R2')\n",
    "plt.xticks(rotation= 90)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Посмотрим на результаты метрик алгоритмов с лучшими параметрами"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "clf = RandomForestRegressor(random_state=17)\n",
    "rf_base = Pipeline([('scalar', scalar), ('RandomForestRegressor', clf)])\n",
    "rf_base.set_params(RandomForestRegressor__n_estimators=150, RandomForestRegressor__max_depth=15).fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = lgb_base.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(rf_base.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Неплохо, но по времени он протигрывает lgb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "clf = lgb.LGBMRegressor(random_state=17)\n",
    "lgb_base = Pipeline([('scalar', scalar), ('lgb', clf)])\n",
    "lgb_base.set_params(lgb__learning_rate=0.06, lgb__n_estimators=900, lgb__max_depth=4).fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = lgb_base.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(lgb_base.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Лучший результат и отличное время"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "clf = Lasso(random_state=17)\n",
    "lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])\n",
    "lasso_base.set_params(Lasso__alpha=0.001).fit(X_train, y_train)\n",
    "                     \n",
    "y_pred = lasso_base.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(lasso_base.score(X_test, y_test), \n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Правда, я не очень верю в стабильность регрессии с alpha = 0.001 - скорее всего на отложенной выборке такая регуляризация будет пагубно влиять на результат, если отложенные данные имеют хоть немного отличное статистическое распределение."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ну и финальная проба - стеккинг и блендинг\n",
    "Воспользуемся авторской реализацие стеккинга и блендинга от Александра Дьяконова отсюда:\n",
    "\n",
    "https://github.com/Dyakonov/ml_hacks/blob/master/dj_stacking.ipynb\n",
    "\n",
    "Отличная статья на эту же тему от Александра Дьяконова:\n",
    "\n",
    "https://alexanderdyakonov.wordpress.com/2017/03/10/cтекинг-stacking-и-блендинг-blending/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Класс для стеккинга и блендинга \n",
    "class DjStacking(BaseEstimator, ClassifierMixin):  \n",
    "    \"\"\"Стэкинг моделей scikit-learn\"\"\"\n",
    "\n",
    "    def __init__(self, models, ens_model):\n",
    "        \"\"\"\n",
    "        Инициализация\n",
    "        models - базовые модели для стекинга\n",
    "        ens_model - мета-модель\n",
    "        \"\"\"\n",
    "        self.models = models\n",
    "        self.ens_model = ens_model\n",
    "        self.n = len(models)\n",
    "        self.valid = None\n",
    "        \n",
    "    def fit(self, X, y=None, p=0.25, cv=3, err=0.001, random_state=17):\n",
    "        \"\"\"\n",
    "        Обучение стекинга\n",
    "        p - в каком отношении делить на обучение / тест\n",
    "            если p = 0 - используем всё обучение!\n",
    "        cv  (при p=0) - сколько фолдов использовать\n",
    "        err (при p=0) - величина случайной добавки к метапризнакам\n",
    "        random_state - инициализация генератора\n",
    "            \n",
    "        \"\"\"\n",
    "        if (p > 0): # делим на обучение и тест\n",
    "            # разбиение на обучение моделей и метамодели\n",
    "            train, valid, y_train, y_valid = train_test_split(X, y, test_size=p, random_state=random_state)\n",
    "            \n",
    "            # заполнение матрицы для обучения метамодели\n",
    "            self.valid = np.zeros((valid.shape[0], self.n))\n",
    "            for t, clf in enumerate(self.models):\n",
    "                clf.fit(train, y_train)\n",
    "                self.valid[:, t] = clf.predict(valid)\n",
    "                \n",
    "            # обучение метамодели\n",
    "            self.ens_model.fit(self.valid, y_valid)\n",
    "            \n",
    "        else: # используем всё обучение\n",
    "            \n",
    "            # для регуляризации - берём случайные добавки\n",
    "            self.valid = err*np.random.randn(X.shape[0], self.n)\n",
    "            \n",
    "            for t, clf in enumerate(self.models):\n",
    "                # это oob-ответы алгоритмов\n",
    "                self.valid[:, t] += cross_val_predict(clf, X, y, cv=cv, n_jobs=-1, method='predict')\n",
    "                # но сам алгоритм надо настроить\n",
    "                clf.fit(X, y)\n",
    "            \n",
    "            # обучение метамодели\n",
    "            self.ens_model.fit(self.valid, y)  \n",
    "            \n",
    "\n",
    "        return self\n",
    "    \n",
    "\n",
    "\n",
    "    def predict(self, X, y=None):\n",
    "        \"\"\"\n",
    "        Работа стэкинга\n",
    "        \"\"\"\n",
    "        # заполение матрицы для мета-классификатора\n",
    "        X_meta = np.zeros((X.shape[0], self.n))\n",
    "        \n",
    "        for t, clf in enumerate(self.models):\n",
    "            X_meta[:, t] = clf.predict(X)\n",
    "        \n",
    "        a = self.ens_model.predict(X_meta)\n",
    "        \n",
    "        return (a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Обучу различные модели с разными гиперпараметрами и посмотрю на их отдельные результаты. \n",
    "\n",
    "df_train = work_data.drop(['id','price', 'date'], axis=1)\n",
    "X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)\n",
    "\n",
    "las0 = Lasso(alpha=0.001)\n",
    "las0.fit(X_train, y_train)\n",
    "y_pred = las0.predict(X_test)\n",
    "print(\"las0, R2_score: {}, RMSLE: {}\".format(las0.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))\n",
    "las1 = Lasso(alpha=0.1)\n",
    "las1.fit(X_train, y_train)\n",
    "y_pred = las1.predict(X_test)\n",
    "print(\"las1, R2_score: {}, RMSLE: {}\".format(las1.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))\n",
    "\n",
    "las2 = Lasso(alpha=0.01)\n",
    "las2.fit(X_train, y_train)\n",
    "y_pred = las2.predict(X_test)\n",
    "print(\"las2, R2_score: {}, RMSLE: {}\".format(las2.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))\n",
    "\n",
    "rf1 = RandomForestRegressor(n_estimators=50)\n",
    "rf1.fit(X_train, y_train)\n",
    "y_pred = rf1.predict(X_test)\n",
    "print(\"rf1, R2_score: {}, RMSLE: {}\".format(rf1.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))\n",
    "rf2 = RandomForestRegressor(n_estimators=250, max_depth=20)\n",
    "rf2.fit(X_train, y_train)\n",
    "y_pred = rf2.predict(X_test)\n",
    "print(\"rf2, R2_score: {}, RMSLE: {}\".format(rf2.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))\n",
    "\n",
    "gbm1 = lgb.LGBMRegressor(boosting_type='gbdt', learning_rate=0.06, max_depth=4, n_estimators=900, nthread=-1, objective='regression')    \n",
    "gbm1.fit(X_train, y_train)\n",
    "y_pred = gbm1.predict(X_test)\n",
    "print(\"gbm1, R2_score: {}, RMSLE: {}\".format(gbm1.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))\n",
    "\n",
    "gbm2 = lgb.LGBMRegressor(boosting_type='gbdt', learning_rate=0.05, max_depth=3, n_estimators=800, nthread=-1, objective='regression')    \n",
    "gbm2.fit(X_train, y_train)\n",
    "y_pred = gbm2.predict(X_test)\n",
    "print(\"gbm2, R2_score: {}, RMSLE: {}\".format(gbm2.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))\n",
    "\n",
    "gbm3 = lgb.LGBMRegressor(boosting_type='gbdt', learning_rate=0.05, max_depth=5, n_estimators=1000, nthread=-1, objective='regression')    \n",
    "gbm3.fit(X_train, y_train)\n",
    "y_pred = gbm3.predict(X_test)\n",
    "print(\"gbm3, R2_score: {}, RMSLE: {}\".format(gbm3.score(X_test, y_test),\n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Набор моделей, в качестве мета-алгоритма возьму Lasso регрессию.\n",
    "models = [las2,las0, rf1, rf2, gbm1, gbm2, gbm3]\n",
    "ens_model = Lasso(alpha = 0.1, random_state = 17)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Запустим блендинг и стеккинг\n",
    "\n",
    "s1 = DjStacking(models, ens_model)\n",
    "s1.fit(X_train, y_train)\n",
    "y_pred1=s1.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(r2_score(y_test, y_pred1), \n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred1,  lmbda_price))))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s2 = DjStacking(models, ens_model)\n",
    "s2.fit(X_train, y_train, p=-1)\n",
    "y_pred2 = s2.predict(X_test)\n",
    "print(\"R2_score: {}, RMSLE: {}\".format(r2_score(y_test, y_pred2), \n",
    "                                       rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(y_pred2,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Стеккинг и блендинг не позволили улучшить качество модели. И в этом нет ничего не обычного, посокльку у нас среди моделей есть очень хорошая модель, которая дает очень высокий результат. Скорее всего улучшений можно было бы добиться больше \"лучших\" разннобразных моделей, которые в итоге бы дали небольшой прирост качства. Но, в целом, практическая польза от такого метода возможна, разве что, в соревнованиях, поскольку с ростом числа объектов в выборке и роста числа базовых тяжелых алгоритмов скорость работы будет еще больше увеличиваться относительно 1 алгоритма, например того же lgb. Еще, как вариант, можно исключать из стеккинга \"тяжелые алгоритмы\" и пробовать добиться высокого качества на простых, но польза от такого опять же только в соревнованиях. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Часть 8. Прогноз для отложенной выборки сравнение оптимальных параметров алгоритма победителя для тестовых и отложенных данных."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Мержим дополнительные таблицы\n",
    "validation_data = validation_data.merge(clinic, how='left', on='zipcode')\n",
    "validation_data = validation_data.merge(shop_horeca, how='left', on ='zipcode')\n",
    "validation_data = validation_data.merge(school, how='left', on ='zipcode')\n",
    "validation_data = validation_data.merge(crime, how='left', on='zipcode')\n",
    "validation_data = validation_data.fillna(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train = work_data.drop(['id','price', 'date'], axis=1)\n",
    "df_test = validation_data.drop(['id','price', 'date'], axis=1)\n",
    "X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Проверяем лучшую модель. Тренируем на всей work_data и предсказываем для validation_data\n",
    "дополнительно немного переберем параметры, дабы понять насколько мы переобучаемся на тестовых данных,\n",
    "когда настраивали гиперпараметры\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "test_R2, holdout_R2, test_rmsle, holdout_rmsle  = [], [],[], []\n",
    "n_estimator = [100,200,300,500,700,900]\n",
    "\n",
    "for i in n_estimator:\n",
    "\n",
    "    lgb_final = lgb.LGBMRegressor(max_depth=4, lgb__learning_rate=0.06, n_estimators=i)   \n",
    "    \n",
    "    lgb_final.fit(X_train, y_train)\n",
    "    test_R2.append(r2_score(y_test, lgb_final.predict(X_test)))\n",
    "    test_rmsle.append(rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(lgb_final.predict(X_test),  lmbda_price)))    \n",
    "    holdout_R2.append(r2_score(target_val, lgb_final.predict(df_test)))\n",
    "    holdout_rmsle.append(rmsle(inv_boxcox(target_val,   lmbda_price), inv_boxcox(lgb_final.predict(df_test),  lmbda_price)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 15))    \n",
    "plt.plot(n_estimator, test_R2, label='test_R2')\n",
    "plt.plot(n_estimator, holdout_R2, label='holdout_R2')\n",
    "plt.title('Ошибка lgb на отложенной выборке R2')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 15)) \n",
    "plt.plot(n_estimator, test_rmsle, label='test_rmsle')\n",
    "plt.plot(n_estimator, holdout_rmsle, label='holdout_rmsle')\n",
    "plt.title('Ошибка lgb на отложенной выборке RSMLE')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Теперь посмотрим на качество при изменении глубины\n",
    "\n",
    "test_R2, holdout_R2, test_rmsle, holdout_rmsle  = [], [],[], []\n",
    "n_estimator = [1,2,3,4,5,6]\n",
    "\n",
    "for i in n_estimator:\n",
    "\n",
    "    lgb_final = lgb.LGBMRegressor(max_depth=i, lgb__learning_rate=0.06, n_estimators=500)   \n",
    "    \n",
    "    lgb_final.fit(X_train, y_train)\n",
    "    test_R2.append(r2_score(y_test, lgb_final.predict(X_test)))\n",
    "    test_rmsle.append(rmsle(inv_boxcox(y_test,   lmbda_price), inv_boxcox(lgb_final.predict(X_test),  lmbda_price)))    \n",
    "    holdout_R2.append(r2_score(target_val, lgb_final.predict(df_test)))\n",
    "    holdout_rmsle.append(rmsle(inv_boxcox(target_val,   lmbda_price), inv_boxcox(lgb_final.predict(df_test),  lmbda_price)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 15))    \n",
    "plt.plot(n_estimator, test_R2, label='test_R2')\n",
    "plt.plot(n_estimator, holdout_R2, label='holdout_R2')\n",
    "plt.title('Ошибка lgb на отложенной выборке R2')\n",
    "plt.legend();\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 15)) \n",
    "plt.plot(n_estimator, test_rmsle, label='test_rmsle')\n",
    "plt.plot(n_estimator, holdout_rmsle, label='holdout_rmsle')\n",
    "plt.title('Ошибка lgb на отложенной выборке RSMLE')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gbm_win = lgb.LGBMRegressor(learning_rate=0.06, max_depth=7, n_estimators=500, random_state=17)    \n",
    "gbm_win.fit(X_train, y_train)\n",
    "y_pred = gbm_win.predict(df_test)\n",
    "print(\"gbm_win, R2_score: {}, RMSLE: {}\".format(gbm_win.score(df_test, target_val),\n",
    "                                       rmsle(inv_boxcox(target_val,   lmbda_price), inv_boxcox(y_pred,  lmbda_price))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Часть 9. Итоговый вывод."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Вцелом за счет преобразования признаков и таргета, а также за счет добавления новых признаков, плюс за счет подбора гиперпараметров модели мне удалось улучшить коэффициент детерминации R2 на 8% и уменьшить среднеквадратичную логарифмическую ошибку почти на 0.06 единиц на тестовой выборке. \n",
    "2. Лучшим алгоритмом обучения по параметру качество - скорость - стабильность стал lightgbm\n",
    "3.  Для lgb качество на валидации немного ниже чем на тесте, но тоже достаточно высокое. Регрессия сильно переобучается и очень не стабиль уже на тестовой валидации. \n",
    "4. Лучшее качество на тестовой выбрке (R2_score: 0.9347084739508015, RMSLE: 0.13765995334788764) с параметрами lr=0.06, md=4, n_es=900, на отложенной(R2_score: 0.9168339280480744, RMSLE: 0.14923617757619653) с параметрами lr=0.06, md=7, n_es=500. Расходжение очень небольшое, что говорит, что lgb очень хорошо себя показал в плане обобщающей способности. \n",
    "5. Точки роста качества мне видятся в поиске новых признаков для данного региона по каждой возможной локации, таких как: средний доход хозяйства, процент жилья сдаваемого в аренду, средняя стоимость проживания, удаленность от административных центоров основных городов, наличие экологических и парковых зон рядом с объектами недвижимости. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
